#Replication Materials, using R 3.61
## Rstan updates frequently, and may have difficulty installing
# For this project I needed to install xcode using  
#terminal to fix bugs in Rstan
rm(list = ls())

#Comment once installed
# install.packages("StanHeaders", dep=T, repos="https://cloud.r-project.org")
# install.packages("rstan", dep=T, repos="https://cloud.r-project.org")
# install.packages("dgo", dependencies=T,repos='http://cran.rstudio.com/')
# install.packages("shape",dependencies=T,repos='http://cran.rstudio.com/')
# install.packages("devtools", dependencies=T,repos='http://cran.rstudio.com/')
# install.packages("dplyr", dependencies=T,repos='http://cran.rstudio.com/')
# install.packages("plyr", dependencies=T,repos='http://cran.rstudio.com/')


library(dplyr)
library(plyr)
library(dgo)
library(dgodata)
library(shape)
library(ggplot2)
library(devtools)
library(StanHeaders)
library(rstan)

setwd("~/Dropbox/Elite Opinion and Labor Policy/Replication Materials")

data.ipoll<-read.csv("data-17may28-nodc.csv")

data.9214<-subset(data.ipoll,YR>="1992")
data.9214<- data.9214[order(data.9214$Abb),]

#Basic Recoding
data.9214$Inc3[data.9214$Income==0|data.9214$Income==1|data.9214$Income==2]<-1 #Under 40k
data.9214$Inc3[data.9214$Income==4|data.9214$Income==5]<-2 #Under 75k
data.9214$Inc3[data.9214$Income==7|data.9214$Income==8]<-3 #Over75

table(data.9214$Inc3)
table(data.9214$Income)

##########THIS IS WHAT WE USED!!!##########
#dgmrp that uses income as a category
data.9214$Year<-data.9214$YR

data.9214$White<-as.factor(data.9214$White)

data.9214$Union.B<-c()
data.9214$Union.B[data.9214$UnionSup>=.51]<-1
data.9214$Union.B[data.9214$UnionSup<=.49]<-0
data.9214$Union.B[data.9214$UnionSup==.5]<-NA

data.9214$Inc3<-as.factor(data.9214$Inc3)
data.9214$Abb<-as.character(data.9214$Abb)


dgirt_in_union <- shape(data.9214, item_names = "Union.B", 
                        time_name = "Year",
                        geo_name = "Abb", 
                        group_names = c("White","Inc3"))

summary(dgirt_in_union)
get_n(dgirt_in_union, by = "Abb")

dgmrp_out_union_2 <- dgmrp(dgirt_in_union, iter = 2500, chains = 4, cores =
                             4, seed = 42)
summary(dgmrp_out_union_2)
head(as.data.frame(dgmrp_out_union_2))
plot_rhats(dgmrp_out_union_2)


#Load data proportions from the CPS
a<-read.csv("popfreq_allyears.csv")
table(data.9214$YR)


a$Abb<-as.character(a$Abb)
a$Inc3<-as.character(a$Inc3)
a$White<-as.character(a$White)

head(a)

b<-aggregate(a$total, by=list(a$Abb,a$YR), FUN=sum)
table(b$x)


dgmrp_out_union_3<-as.data.frame(dgmrp_out_union_2)
dgmrp_out_union_3$YR<-dgmrp_out_union_3$Year

ps_income <- poststratify(dgmrp_out_union_3, a, strata_names =
                            c("Abb", "YR","Inc3"),proportion_name="total",
                          aggregated_names = c("White"))

head(ps_income)
write.csv(ps_income, file = "dgmrp-union-raceinc-RR.csv")


d<-dgirt_plot(ps_income, group_name = "Inc3", time_name = "YR", geo_name = "Abb")
head(d)
d$data
g<-d$data
write.csv(d$data,"Median-CIs-States.csv")

theme_set(theme_bw())

#Black and white plots
pdf(file="sub_dgmrp.pdf",width=8,height=7)
plot<-ggplot(d$data , aes(x = YR , y = median, fill=Inc3)) +
  facet_wrap(~Abb) +
  geom_ribbon(data=d$data, aes(ymin=q_025, ymax=q_975, group = Inc3), alpha = 0.35)+
  theme(axis.text.x=element_text(angle=90))
plot<-plot+ scale_x_continuous(name="Year", limits=c(1992, 2014))
plot<-plot+ ylab("Union Support")
plot<-plot + scale_fill_manual(values=c("#CCCCCC", "#666666", "#000000"),
                               name="Income Third",
                                 breaks=c("1", "2", "3"),
                                 labels=c("Lowest", "Middle", "Highest"))
plot
dev.off()

# Black and White Subplot
 h<-subset(g,Abb=="OH"|Abb=="SC"|Abb=="CA"|Abb=="NY")

pdf(file="sub_dgmrp_1.pdf",width=8,height=5)
sub<-ggplot(h , aes(x = YR , y = median, fill=Inc3)) +
  facet_wrap(~Abb) +
  geom_ribbon(data=h, aes(ymin=q_025, ymax=q_975, group = Inc3), alpha = 0.35)+
  theme(axis.text.x=element_text(angle=90))
sub<-sub+ scale_x_continuous(name="Year", limits=c(1992, 2014))
sub<-sub+ ylab("Union Support")
sub<-sub + scale_fill_manual(values=c("#CCCCCC", "#666666", "#000000"),
                               name="Income Third",
                               breaks=c("1", "2", "3"),
                               labels=c("Lowest", "Middle", "Highest"))
sub
dev.off()



# #Color Plots
# pdf(file="sub_dgmrp_color.pdf",width=8,height=7)
# plot<-ggplot(d$data , aes(x = YR , y = median, fill=Inc3)) +
#   facet_wrap(~Abb) +
#   geom_ribbon(data=d$data, aes(ymin=q_025, ymax=q_975, group = Inc3), alpha = 0.3)+
#   theme(axis.text.x=element_text(angle=90))
# plot<-plot+ scale_x_continuous(name="Year", limits=c(1992, 2014))
# plot<-plot+ ylab("Union Support")
# plot<-plot + scale_fill_discrete(name="Income Third",
#                                  breaks=c("1", "2", "3"),
#                                  labels=c("Lowest", "Middle", "Highest"))
# plot

# #Color Subplot
# h<-subset(g,Abb=="OH"|Abb=="SC"|Abb=="CA"|Abb=="NY")
# 
# pdf(file="sub_dgmrp_1_color.pdf",width=8,height=5)
# sub<-ggplot(h , aes(x = YR , y = median, fill=Inc3)) +
#   facet_wrap(~Abb) +
#   geom_ribbon(data=h, aes(ymin=q_025, ymax=q_975, group = Inc3), alpha = 0.3)+
#   theme(axis.text.x=element_text(angle=90))
# sub<-sub+ scale_x_continuous(name="Year", limits=c(1992, 2014))
# sub<-sub+ ylab("Union Support")
# sub<-sub + scale_fill_discrete(name="Income Third",
#                                breaks=c("1", "2", "3"),
#                                labels=c("Lowest", "Middle", "Highest"))
# sub
# dev.off()





